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Electronic energies are calculated for a Hubbard model on the Ceo molecule using projector 
quantum Monte Carlo (QMC). Calculations are performed to accuracy high enough to determine 
the pair binding energy for two electrons added to neutral Ceo- The method itself is checked against a 
variety of other quantum Monte Carlo methods as well as exact diagonalization for smaller molecules. 
Fhe conclusion is that the ground state with two extra electrons on one Ceo molecule is a triplet 
and, over the range of parameters where QMC is reliable, has a slightly higher energy than the state 
with electrons on two separate molecules, so that the pair is unbound. 

PACS numbers: 74.70.Wz, 71.10.Li, 02.70.Ss 



I. INTRODUCTION 

The discovery of superconductivity in the alkali- metal- 
doped bulk fullerenes K 3 Ceo and Rb^C^^r^ sparked in- 
tense interest in fullerene-based materials, leading to ex- 
tensive experimental and theoretical studies (for a review 
see Ref. Q) ■ Theoretical calculations to explain the insu- 
lating and superconducting properties of bulk fullerene 
materials fall into two major categories: molecular level 
calculations which determine the effective interactions 
of intramolecular elect rons^£*§*£*£*2*i2iii and lattice level 
calculations based on an effective Hamiltonian in which 
the intramolecular degrees of freedom have been inte- 
grated out. 12 i 13 i 14 

Although much of the work has focused on the phonon 
mechanism for superconductivity in the alkali-Cgo's, 
some authors have proposed a purely electronic mech- 
anism. In particular, it was argued by Chakravarty, 
Kivelson and co-workers (CK) that electronic interac- 
tions within a single 6*60 molecule can lead to an effective 
attraction between charge carriers This argu- 

ment was supported by perturbative calculations of the 
electronic binding energies of the conventional one band 
Hubbard model, on the Cgo structure. Results of the 
CK calculation suggest that electrons have a tendency to 
form paired states in a single fullerene molecule, rather 
than remain separate. This tendency could be the origin 
of the attractive interaction which is an essential ingre- 
dient of the BCS theory of superconductivity. 

However, one might doubt the applicability of pertur- 
bation theory to this problem. First, the Hubbard repul- 
sion U in the CK calculation is approximately 75% of the 
bandwidth, so it is hardly a small parameter. Also, the 
binding energy is typically a small quantity, calculated 
from the difference of large internal energies of the Cqq 
molecule at different electron dopings. Low order per- 
turbation theory estimates of such subtle energy differ- 
ences may be unreliable. Thus it is interesting to repeat 
the calculation using different methods, which might lend 
support to or cast doubt upon the perturbation theory 
results. 

In this paper we use quantum Monte Carlo (QMC) 
calculations to estimate the binding energy of pairs of 
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FIG. 1: Comparison of electronic pair binding energies 
Ab(61)/t, defined in Eq. (|5J, obtained from perturbation the- 
ory in different spin sectors (solid and dash-dot lines)-*"!— and 
PQMC calculations on a Ceo molecule. PQMC finds S = 1 
for ground states with 62 electrons. 



electrons on a single Ceo Hubbard molecule. In order to 
establish a high level of confidence in our results, we use 
a number of complementary QMC methods, including 
auxiliary field QMC for both real and imaginary chemi- 
cal potentialsi^iiS and stochastic series expansion (SSE) 
at finite temperature T, and projector QMC (PQMC) at 
T = on a series of Hubbard molecules, with the number 
of sites ranging from 4 to 60. In addition, comparison is 
made to the results of exact diagonalization (ED) for sys- 
tems of up to 12 sites. Our main result, shown in Figure 
n is a comparison of the pair binding energy (see below 
for a definition) for two electrons added to a neutral Cqq 
molecule, to the perturbation calculations of CK. All en- 
ergies are measured in units of the hopping parameter 
t. In contrast to perturbation theory, which finds that 
the ground state is a singlet for U/t > 3, we find, in 
agreement with Hund's rule, that the ground state re- 
mains a triplet state over the entire range of U studied. 
In particular, there is no indication of the attractive sin- 
glet ground state which perturbation theory finds for U/t 
greater than about 3.3. Our QMC studies find small pos- 
itive binding energies for U/t < 4.5, indicating that two 
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TABLE I: Comparison of exact diagonalization and PQMC 
calculations on the truncated tetrahedron (12 sites) at U = 2t. 
PQMC simulation parameters: /3 = 10/t, At = 0.05/t, N m = 
10 7 . E n (Sz) is the energy of a system with n electrons and 
z-component of total spin S z . A„, m is the energy difference 
Ei2+ n {S") - E 12+ m(ST) with {st,ST) given in the second 
column. For binding energies Ab(n) the second column shows 
(S z +1 , 5*" -1 , S z ) - the S z values for 3 states involved in its 
calculation, in the order of appearance in Eq. (2) 



separate molecules, each with one extra electron, have 
lower energy than one molecule with two extra electrons. 
The largest value of U we are able to study is 4.5i, since 
the sign problem discussed below becomes unmanageable 
for larger values of U. 

The paper is organized as follows: in the next section 
we introduce the model and the QMC methods used in 
our simulations. We then proceed to present the tests of 
QMC codes on smaller molecules (truncated tetrahedron, 
etc.), where exact analytical or exact diagonalization re- 
sults are available. Then the results for the Cq$ molecule 
are presented and analyzed. Finally, a conclusion based 
on our numerical results is drawn and the reliability of 
the method is discussed. 



II. METHODOLOGY 

Following CK, we consider a one-band Hubbard model 
with the Hamiltonian H = Hq + Hi, defined on a Cqq 
molecule by 

Hq = ~ y~] tij(cl a c ja + h.c.) - ^y^nio-, 

(ij)tr i<r 

Hi = U^n^riii - — ^n iCT . (1) 

i ia 

Here H contains the standard kinetic energy and chemi- 
cal potential term. The summation in the kinetic energy 





U = 


n 


-1 





1 


^exact 


0.7959095 


1.0001765 


1.2043146 


rp 

-L^ exact 


-44.2672020 


-46.3708440 


-44.3478000 


"•AFQMC 


0.7959095 


1.0001765 


1.2043146 


-tiAFQMC 


-44.2672020 


-46.3708440 


-44.3478000 


OAFQMC 


1.0 


1.0 


1.0 




U = At 


n 


-1 





1 


"SSB 


0.873(1) 


1.0005(2) 


1.126(1) 


Esse 


-16.5(2) 


-4.0(2) 


13.9(1) 


Ssse 


0.955 


0.957 


0.960 


flAFQMC 


0.8734(2) 


1.000078(1) 


1.1266(2) 


-Eafqmc 


-16.61(2) 


-4.13(4) 


13.74(6) 


Safqmc 


1.0 


1.0 


1.0 



TABLE II: Comparison of the density n, total internal en- 
ergy E and average sign S between exact analytical results 
(U = 0), SSE (U = At) and AFQMC on a C* 60 molecule. Sim- 
ulation parameters: f3 = 0.5/t, At = 0.05/i, N m = 10 s . In 
the SSE run 10 7 measurements, separated by a full diagonal 
and directed loop update, were performed. At lower tem- 
peratures SSE is unreliable due to the severe sign problem. 

term is performed over all nearest neighbor pairs on a 
Cqq molecule. The hopping constants Uj are chosen to 
be equal to t for the single bonds, connecting a pentagon 
and a hexagon, and to t' = 1.2t for the double bonds be- 
tween two hexagons. Hi is a sum of the on-site Coulomb 
repulsion (Hubbard) term and a diagonal term, added to 
make the model particle-hole symmetric around /i = on 
bipartite lattices. Clearly, this additional term does not 
affect the value of the electronic binding energy, which 
we choose to define as 

A b (n) = £ n+1 + £ n _ !-2£„, (2) 

where E n is the internal energy of a molecule with n 
electrons. Note that this definition has the opposite sign, 
compared to that of CKi4 In our case the tendency of 
the electrons to bind into pairs is indicated by a negative 
value of the binding energy At,. 

Determinant or auxiliary field QMC (AFQMC) has 
been widely used in model Hamiltonian simulations since 
its introduction by Blankenbecler et al.M**± and its 
further development by Hirsch^SiS^ and White et alm& 
The application of this technique to the one-band Hub- 
bard model starts with the Suzuki- Trotter discretization 
of the imaginary time in the grand canonical partition 
function^ 

= Tre _/3(Ho+Hl) 

= Tvf[e- A ^ H °+ H ^, (3) 

8=1 

where (3 — l/(fcsT) is the inverse temperature, dis- 
cretized in such a way that — AtL. After the applica- 
tion of the Hubbard-Stratonovich transformationS&SL the 
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FIG. 2: Grand canonical simulation of a C%o molecule at various chemical potentials. Simulation parameters: U = At, f3 = 10/t, 
At = 0.1/t, Nt — 10 3 , N m = 10 5 . We are only interested in the qualitative behavior of the system around half-filling, so the 
statistical errors are not estimated, (a) Energy per site vs. chemical potential, (b) Electron density vs. chemical potential, 
(c) Average sign as a function of chemical potential, (d) Average sign as a function of electron density. Curves connecting the 
points are guides to the eye. 



fcrmionic degrees of freedom in Eq. may be traced 
out, and we arrive at 

II det[l + B L (a)B L ^(a) ■ ■ ■ B 1 (a)} 
= ^detO({<r},/i) T detO({<7},/i) i . (4) 
The Bi matrices are defined as 

Bl {a) = e -±rK/2 e V°W e -ArK/2 } (5) 

f —Uj for i,j nearest neighbours, 
[K)z: > ~ \ otherwise, (bj 

V$(l) = %[Ac^(0 + mAt], (7) 

where (T,(Z) = ±1 is the auxiliary Ising spin coupled with 
the electrons at lattice site i and time /At, and a = ±1 
corresponds to f or J, in Eq. J3J|. In Eq. JSJ we have 
used a symmetric decomposition of the partition func- 
tion, which produces a much smaller Trotter erro»22i2i 
compared to the non-symmetric decompositioni22i2&2!i2& 
The Monte Carlo (MC) weight P is then given by the 
product of two determinants in Eq. ( which is always 
positive for bipartite lattices at half-filling>S At low tem- 
peratures we use QR factorization to stabilize the matrix 



multiplications and inversions'^ A version (projector 
QMC or PQMC) of the above procedure can be used to 
directly project out the ground state properties from an 
initial trial wave function; see Ref. |2fj for details. 

Unless indicated otherwise, PQMC and AFQMC cal- 
culations were performed at the temperature fixed by /3 = 
10/t, with imaginary time discretization At = 0.05/i. 
The system was first brought to thermal equilibrium by 
performing Nt thermalizations sweeps S2, followed by 
N m measurements with a single sweep SI— performed 
between them. 

Since we are interested in non-bipartite molecules such 
as Ceo, the MC weights in general are not always positive. 
In the case of negative weight P, we associate a proba- 
bility value \P\ with it, and include a sign S = P/\P\ in 
the average: (E) = (ES)/(S). Now the average (. . .) is 
with respect to the probability distribution \P\. 

Finally, we note that in estimating the statistical error 
for a composite quantity X = X\ + X 2 + ■ ■ ■ + X n , we use 
the standard formula SX = (SXf + 8X% + ... + SXl) 1 / 2 , 
where SXi's are estimated by the Jackknife method— 
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TABLE III: PQMC calculations on a 6*60 molecule. Part A of the table shows the total internal energy E n (S z ) of a Ceo 
molecule with n electrons and z-component of total spin S z . The parameters used in the simulations are t' = 1.2t, f3 = 10/t, 
At = 0.0625/t (for U = 4t), At = 0.05/t (for other U values), N m = 10 7 . N m data were divided into 10 bins for error 
estimation. For n — 60, 61, 62, we have collected more data (between 4 x 10 7 and 8 x 10 7 measurements) for a more accurate 
comparison between PQMC and perturbative results. Part B shows electron (hole) binding energies Ab(n). As before, the 
S z column in this case lists the S z values of three states, involved in the calculation of the binding energy, in the order of 
appearance in Eq. Q. For example, A b (58) with S z = (1/2,3/2, 1) denotes E 59 (S Z = 1/2) + £ , 57 (3/2) - 2E 5 $(1). The data 
points marked with * were calculated using a non-symmetric decomposition in Eq. Jjjjt. Only limited results were obtained for 
U — 4.5t because of the long averaging times required. 



III. APPLICATION 
A. Comparison to other methods 

In this section, we check our QMC programs against 
ED and SSE results. 34 i 35 i 36 Table [Q lists energies and 
binding energies from both ED and PQMC for a trun- 
cated tetrahedron, which has 12 lattice sites and 3 nearest 
neighbors for each site. For energies at different dopings, 
good agreement is obtained between the two methods. 
The largest deviation of PQMC from ED is found for 
n = 10 (S z = 0) (about 0.7% deviation), which might be 
due to the relatively low value of the average sign and 
the incomplete projection of a nearby singlet excitation. 
Energy differences between two different dopings (e.g., 
Ai : o) are in good agreement for the two methods. How- 
ever, the inaccuracies are magnified when pair-binding 



energies are extracted from two already-small energy dif- 
ferences, although some of these energies are still in good 
agreement for the two methods within error bounds, e.g., 
A b (ll). The difficulty of extracting A b (13) from PQMC 
is possibly because the ground state with 14 electrons lies 
in the spin singlet sector (as confirmed by ED), and it is 
difficult for PQMC to completely project out the nearby 
spin triplet state (first excited state). There is no such 
problem if the ground state for two-electron doping is a 
spin triplet, which, as we will see below, is exactly the 
case for Cqq . From the good agreement between ED and 
PQMC, we conclude that the discretization error caused 
by At = 0.05 /t is sufficiently small. We also find that 
the projection factor f3 = 10/t is large enough to project 
out the ground state from an initial trial state. We will 
use these values of Ar and (3 in our AFQMC and PQMC 
simulations of the Cgo molecule. 
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In Table El we check our grand canonical simulation 
program (AFQMC), against ED (at [7 = 0) and SSE (at 
U = At) on a 6*60 molecule. Again we see good agree- 
ment in both the density n and energy E calculations 
among these methods. For {7 = 0, the AFQMC results 
are exactly the same as the exact diagonalization results. 
This is because at U = 0, there is no coupling between 
electrons and the auxiliary Ising field; the Ising field is 
wiped out completely and the electrons cannot feel the 
existence of Ising spins. The simulation at U — also 
shows that the discretization error is absent in AFQMC. 
In both U = and U = At, we have set the temperature 
T = 2t to avoid a severe sign problem in the SSE simula- 
tion. Because SSE does not suffer from the discretization 
error, we again confirm that At = 0.05/t in AFQMC is 
sufficiently small to avoid any systematic discrepancy. 

The results in Table U for simulations on a truncated 
tetrahedron molecule show that PQMC results are in 
good agreement with ED. The systematic discretization 
error caused by At is reasonably small; thus the extrap- 
olation to At = is unnecessary. 



B. Application to Ceo molecule 

In this section, we discuss the QMC results for a Cqo 
molecule. Figure shows the results of an AFQMC sim- 
ulation of the Hubbard model on a Ceo molecule at var- 
ious chemical potentials \i. At half-filling, unlike the 
bipartite 2D square lattice, the QMC simulation has a 
slight sign problem due to the pentagon frustration in 
the Cgo geometry; see panels (c) and (d) in Figure |21 
From (c) and (d) we also see that hole dopings have a 
worse sign problem than electron dopings. As expected, 
the compressibility^ k = ^r^ 2 ~ around /i = 0. 

Tabic EH lists the PQMC results for 2<U/t< 4.5 for 
the Cqo molecule. It can be seen that hole doping causes a 
more severe sign problem than the electron doping, which 
is consistent with the behavior in Figure [21 In Part A 
of the table we see that reasonably accurate results can 
be obtained for U — 2t, 3t and At. The sign problem 
quickly becomes severe for U > At, as is evident in the 
data presented for U — 4.5t, where the average sign is 
only 0.18 for 62 electrons. For U — 5t and 62 electrons, 
the average sign is 0.08 and it is not possible to extract 
a reliable binding energy. 

From TablclTTTR we can calculate the energy difference 
between different total spin sectors at the same doping to 
compare the ground state spin configuration with Hund's 
rule. Let us first discuss the non-interacting single elec- 
tron energy levels of a Cqo molecule. At half-filling, 60 
electrons move independently in "tt" molecular orbitals 
formed from the 60 p z atomic orbitals of the 60 carbon 
atoms. An exact diagonalization of the non-interacting 
Hamiltonian gives 60 energy levels, of which the lowest 30 
levels are occupied at half filling. The lowest unoccupied 
molecular orbitals (LUMO) are 3-fold degenerate. The 
highest occupied molecular orbitals (HOMO) are 5-fold 
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FIG. 3: Huckel energy level diagram for the neutral C§o 
molecule. The lowest 30 levels are doubly occupied. The 
energy level scale is drawn according to the exact diagonal- 
ization of the non-interacting on-site Hubbard Hamiltonian, 
i.e., U — 0. The energy level labels are from those of the icosa- 
hedral group, which is the symmetry group of a C^o molecule. 
The LUMO band is labelled by t lu , and HOMO by h u . We 
will consider doping of LUMO and HOMO for a discussion of 
Hund's rule. 



degenerate. The detailed energy levels of the neutral Cqq 
molecule are shown in Figure|31 The energy gap between 
HOMO and LUMO is l.OAt from this exact diagonaliza- 
tion for non-interacting neutral Cqq- 

We calculate the energy difference AE of different to- 
tal S z values at the same filling as shown in Figure 0| As 
a severe sign problem occurs for U = 5i, the statistical 
errors are too large for a definite discussion of Hund's 
rule. Thus, we mainly discuss results in Figure^] (a)- (c). 
At n — 59,60 and 61, AE is positive, which can be ex- 
plained by the non-interacting energy levels in Figure 
For example, for n — 59 and total spin S z = 1/2, we 
need to flip one electron from spin down to spin up in or- 
der to get S z = 3/2, which means we need to excite one 
spin-down electron from the HOMO band to the LUMO 
band, with an energy cost of 1.04i. The same explana- 
tion applies to n = 60 and 61. At n = 57, 58, 62, 63, 64, 
AE is close to zero which means that the two different 
values of S z are part of the same multiplet. For these 
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FIG. 4: Comparison of the PQMC spin configuration of a Ceo 
molecule at various dopings with Hund's rule. AE is defined 
as the energy difference between two total spin z component 
sectors, e.g., A£(60) = E 60 (S Z = 1) - E 60 (S Z = 0) for the 
neutral molecule and A£(61) = E 61 (S Z = 3/2) - E 61 (S Z = 
1/2) for one-electron doping. A positive AE at fillings n — 
59, 60, 61 can be understood in the non-interacting picture in 
Figure |3] and a negative AE at fillings n = 57,58,62,63,64 
are in agreement with Hund's rule, n = 65, 66 can again be 
explained with Figure [3] See text for discussions. The dotted 
lines connecting MC points are only guides to the eye. 
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state is a triplet, as for n = 62, and if the first excited 
state is a singlet, then the S z = state may contain a 
small admixture of the excited singlet state, while the 
S z = 1 state will not. This will result in the S z = state 
having a slightly higher energy than the S z = 1 state, 
providing a measure of the effectiveness of the projection 
within this multiplet. 

At n = 65, 66, the positive AE can again be explained 
by the single electron picture. A spin-down electron in 
the t lu band must flip its spin and be excited to the 
tig band, with an energy cost of 0.2603t for the non- 
interacting Hamiltonian. This is consistent with Figure 
0]for n = 65 and 66. 

In Table HUB we calculate the pair binding energy 
Ab(n) at various dopings. In these calculation, we have 
used the lowest energy of the different total spin z states 
for a given n. Figure ^ shows a comparison between 
the perturbation calculations of Refs. |j|5| and PQMC 
calculations. There is no indication from the PQMC cal- 
culations of a bound singlet state for U > 3.3i as sug- 
gested by perturbation theory, represented by the solid 
line in Figure^] Instead, the ground state for 62 electrons 
is a triplet over the entire range of parameters studied, 
in agreement with Hund's rule. Furthermore, we find 
that this triplet state is unbound for U/t < 4.5. Un- 
fortunately, the sign problem precludes us from studying 
larger values of U using PQMC. We have also checked the 
energy difference obtained from PQMC with the results 
from imaginary chemical potential simulations, which 
shows good agreement and will be presented elsewhereJ^ 



IV. CONCLUSIONS 

We have performed extensive QMC simulations on a 
single Cqq molecule. The PQMC simulation calculates in- 
ternal energies at various fillings, and shows that Hund's 
rule is well obeyed. In contradiction to the perturbation 
theory result^*^^ we find no singlet pair binding (i.e., 
no negative Ab(n)) for the parameter ranges explored 
(U = 2t,3£,3.5£,4i,4.5t, *' = !- 2 *)- Therefore, a purely 
electronic attractive interaction, originating from the one 
band Hubbard model with on-site Coulomb interaction, 
seems unlikely. This main result is presented in Figure 
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cases, in agreement with Hund's rule, the electrons tend 
to occupy the degenerate HOMO or LUMO in a way that 
maximizes the total spin. 

Note that the small differences in energy for fixed n 
and different S z within a multiplet must result from small 
admixtures of excited states. For example, if the ground 
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